


oo 



cd 



O 
o 



c^ 



Interaction of a granular stream with an obstacle 

Volkhard Buchholtz and Thorsten Poschel* 
Humboldt- Universitdt zu Berlin, Institut fiir Physik, 



D . Invalidenstrafie 110, D-10115 Berlin, Germany 



(February 1, 2008) 



(p ■ Abstract 



C^ . We investigate numerically the interaction of a stream of granular particles 

with a resting obstacle in two dimensions. For the case of high stream velocity 

Vh ■ we find that the force acting on the obstacle is proportional to the square of the 

^' 

P^ . stream velocity, the density and the obstacle size. This behavior is equivalent 

o, 

, ^, . to that of non-interacting hard spheres. For low stream velocity there appears 

CN . a gap in between the obstacle and the incoming stream particles which is filled 

> 
-^ . with granular gas of high temperature and low density. As soon as the gap 

(N 

^ . appears the force does not depend on the square of velocity of the stream but 

C"^^ . the dependency obeys another law. 

0\ 



I. INTRODUCTION 

When a stream of solid frictionless particles i with radii i?j -^ which do not interact 
with each other is scattered at an obstacle one expects from momentum conservation ||l| 



X 

\-{ ' that the net force F acting on the obstacle obeys the law 



F^p-A-vl. (1) 
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A is the projection of the obstacle to the stream direction, Vst is the velocity of the stream 
particles, and p is the density of the stream. When the stream consists of granular particles 
which do interact with each other, i.e. between which forces and torques are transferred, 
and which undergo damping during collisions one can expect deviations from the behavior 
in eq. |I|. The aim of this paper is to investigate the interaction between the stream and the 
obstacle using soft particle molecular dynamics in two dimensions. As far as we know the 
described phenomena have not been investigated experimentally so far. 

The paper is organized in three parts. In section || we explain the numerical experiment 



in detail. Section III A contains our results concerning the angular distributions of the 



scattered stream, i.e. the distributions of mass, momentum and angular momentum as 



functions of the scatter angle G. In IIIB we investigate the force experienced by the obstacle 



from the impacting granular stream, whereas in section |I^ we present some continuum 
properties of the experiment, i.e. the fields of granular temperature and density. In the final 
section |V] we discuss the results. 

II. THE NUMERICAL EXPERIMENT 

The setup of the numerical experiment is sketched in fig. |l|, b is the width of the stream, 
Vst the velocity of the stream particles. When the obstacle is a sphere we can describe it 
by its radius Rob- At distance Ro from the center of the obstacle with Roh <^ Rd there is 
a screen where we measure the angle dependent distributions of the scattered stream. In 
the simulation we used Rd > 12 ■ Rob- The radii of the stream particles Ri are equally 
distributed in the interval [0.05; 0.11] cm. Each of them has translational and rotational 
degrees of freedom. 

Before interacting with the obstacle the particles move in an uniform stream (left in 
fig. |1]), all with the same initial velocity Vst- When released the particles move towards the 
obstacle, collide with the obstacle and once cross the screen plane (right in fig. |1|). After 
crossing the screen plane the particles are removed from the calculation. 



Initially the stream particles are arranged in columns. All particles which belong to the 
same column have the same x-position. One particle of each column has been set to a random 
y-position, then the column was filled in positive and negative y-direction up to ymin = —b/2 
and r/max = b/2 respectively. At the initialization time the distances between the centers of 
two neighboring spheres with radii Ri and Rj in horizontal and vertical direction Ax and 
Ay are 

Ax = 2 ■ R'^^'^^p (2a) 



Ay = {2- R'^"'' + 0.2 * GRND) ^1/p , (2b) 

where p is the density of the stream as used in eq. |l|, R'^°'^ is the radius of the largest 
stream particle and GRND is a Gaussian random number with standard deviation of one. 
All particles move with the same initial stream velocity Vgt in positive x-direction. The 
density value p = 1 in eqs. |^ means that the distance between neighboring columns is 
2 Ji'max^ According to the described method the number of particles is not conserved during 
the simulation. Particles are removed from the calculation when they cross the screen plane, 
and particles are added to the system according to the initialization scheme. The total 
number of particles which are involved in the simulation varies between N = 4000 and 
N = 6000. 

Two particles i and j at positions f^j and fj interact if the distance between their center 
points is smaller than the sum of their radii. In this case we define the particle overlap 
^ij = Ri + Rj — \fi — fj\. Two colliding spheres feel the force 

F,, = F,^.n^ + Ff^.n' (3) 



with the unit vectors in normal direction n^ = Jj Z\ and in shear direction n"^, 



n^ 



1, 



n^ -n^ = 0. 

In viscoelastic approximation the normal force acting between colliding spheres has been 
derived rigorously by Brilliantov et al. [H 



^\/^i, /2 , r -. 

Fl? = 7T^ k^' + AxT^^) (4) 



with 






(5) 



3 (3r/2 + 2r7i) 

y, V and 771/2 are the Young- Modulus, the Poisson ratio and the viscous constants of the 
particle material. R^^--' is the effective radius of the spheres i and j 

,eff_ Ri-Rj 



R 



V 



Ri + Rj 



(6) 



The functional form of eq. (|^) has been found before by Kuwabara and Kono using 
heuristic arguments and has been derived recently by Morgado and Oppenheim [Q from first 
principle calculations. 

For the shear force we apply the Ansatz by Haff and Werner [^j which has been proven 
to give correct results in many other molecular dynamics simulations of granular material, 



e.g. 



sig 



:n (^.70 



mm <; 75 ■ m^f ■ 



rel 



^^■ 



7N 



(7) 



The relative velocity v[ between the particles i and j with angular velocities f2j and VLj 



eff 
and the effective mass rrij^- of the particles i and j are given by 



^' = ffj - fj) ■n'^ + Ri-VLi + Rj ■ Qj , 



771 



ejf. 



TTii ■ rrii 



rrii + rrij 
The resulting torques Mi and Mj acting upon the particles are 



(8) 
(9) 



Mi = F'-Ri , M, 



u 



-^ij ■ ^J 



(10) 



Eq. (0) takes into account that the particles slide upon each other for the case that the 
inequahty 



l^-\F, 



N 



< 



U 



(11) 



holds, otherwise they feel some viscous friction. 



The particles interact with the spherical obstacle in the same manner as the particles 
with each other, i.e. the force Fij acting on the zth particle when colliding with the obstacle 
is given by eq. (^, with Rj = Rob and Vj is the position of the obstacle. The mass of the 
obstacle is infinite, hence the effective mass (eq. (^) is m"^" = rrii. The fiat obstacle was 
built up of equal spheres as shown below in fig. |^ with the same material constants as the 
stream particles but of infinite mass. Their radii have been i?*^ = Ri, where Ri is the mean 
radius of the stream particles. The force acting on a stream particle when colliding with 
such a particle is calculated according to eq. (|). 

For the integration of Newton's equation of motion we used the classical Gear predictor- 



corrector molecular dynamics algorithm |lTO| , pJ]| of 6th order for the positions and of 4th 
order for the rotation of the particles. 

In the numerical simulation we used the parameters given in Table | which are valid 
for a typical granular material. The time step for the numerical integration of the Gear 
predictor-corrector scheme was At = 1 ■ 10~^ sec. Using this time step the system behaves 
numerically stable, i.e. tests with time step At = 2 ■ 10^^ sec lead to the same results. 

The simulations presented in this paper have been performed on the 16 processor message 
passing parallel computer KATJA ( |iittp://summa. physik.hu- berlin.de/KA'l' J A/| ). 



III. RESULTS 

A. Scattering of the granular stream 

First we investigate the angle dependent distributions of the scattered particles for a 
numerical experiment with a fixed spherical obstacle of radius Rob = 1.5 cm. The differential 
mass intensity dm of the scattered stream normalized by the incoming stream is defined by 

dm{e) = ^"^^^Q^ , (12) 

mst 

where dm£,{Q) is the mass of the particles which are scattered into the interval [O, O + dQ] 
and rrist is the total mass of the incoming particles. B = is the direction directly behind 



the obstacle. 

Figure ^ shows the differential intensity, i.e. the mass dm{Q) of the particles which are 
scattered into a certain interval [0,0 + A0], as a function of the scatter angle 0. The 
problem is symmetric due to the symmetry of the incoming particle stream, hence the 
intensity is symmetric as well. For the sake of better statistics we plot in fig. ^ dm{Q) + 
dm{—Q) over 0. We find maxima at ~ ±0.5 = 29°. The minimum at ~ is due to the 
finite size of the obstacle, so only few particles are scattered into the region directly behind 
the obstacle. 

Due to the different sense of the transferred angular moments resulting from collisions 
at the lower or at the upper side (in y-direction) of the particles the angular momentum 
L(0) of the particles which fall into a certain interval of scatter angle + dQ has a negative 
sign for < and a positive for > 0. Hence, the distribution is antisymmetrically. 
As in the case of the intensity we find symmetry of the absolute value of the distribution 
of the angular velocity. Therefore, we draw in fig. ^ -^(0) — L{—Q) over 0. When the 
particles collide with the obstacle they lose kinetic energy due to damping. The larger the 
distance \y\ oi the particles in vertical direction from the obstacle (see fig. |I|) the smaller is 
the dissipated kinetic energy. Hence, the particles moving with large y (outer particles) are 
faster than others with lower \y\ (inner particles). Therefore, the stream particles undergo 
more intensive collisions from the outer side (large \y\) than from the inner (smaller |?/|). 
Particles scattered into the region behind the obstacle have a higher angular velocity since 
they need to be pushed by faster particles multiple times from the same side to reach this 
region. These effects are the reason for the distribution L{<d) of the angular momentum. 
In agreement with the above arguments fig. ^ shows that the distribution of the absolute 
value of the averaged angular momentum L{Q) of particles decreases with increasing 0. In 
the range 0^0 the opposite rotational sense for particles from the upper and from the 
lower side combine to an angular momentum of zero. In this region the error bars are large 
according to the bad statistics. 



(>f')%(4"T 
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The averaged absolute value of the momentum ImjiTjl = m. 
scattered particles as function of the scatter angle is drawn in fig. ^. The velocity of the input 
stream is Vst = 500 cm/sec. As described above the energy loss of the stream is higher for 
particles scattered by low angles G. In this region we find large error bars according to bad 
statistics - only few particles are scattered directly behind the obstacle. The maxima of the 
distribution p(0) are close to the intensity maxima. There seems to be a preferred region for 
the scatter angle where the traces of the particles are almost in parallel, i.e. where energy 
dissipation is low. We will come back to this presumption when we discuss the granular 
temperature field of the stream. 

B. Force experienced by the obstacle 

The impact of the stream particles on the target results in an effective force between the 
stream and the obstacle. With the strong simplifying assumptions that the particles have no 
rotational degree of freedom and the particles interact only with the obstacle but not with 
each other one obtains from kinetic gas theory a theoretical approximation for the force F 

F = f^r^Ro,pvl. (13) 

The force is the transferred momentum per time which is the mass of the particles colliding 
with the obstacle in the time dt times the stream velocity. Eq. (p!3D implies the relations 



F ~ p , (14a) 

Fr^Rob, (14b) 

and F r^vl^ . (14c) 

In the following we will investigate each of the relations ( p!^ numerically for granular 
particles for which the simplifying assumptions are evidently not provided. 

The force acting on the obstacle fluctuates according to the irregularities of the stream, 
hence, to get reliable values we averaged the force for each data point over a certain time 



that corresponds to approximately 500 rows of stream particles colliding with the obstacle. 
As shown in the results this time is sufficient to generate very small error bars. 

As we see from fig. ^ the simulation is in good agreement with eq. ( |14a| ) for the force F 
as a function of the stream density p of the input stream. The function F = F{p) can be 
fitted to 

cr Qjji 

F = (125000 ■ p + 1400) ^^ . (15) 

sec^ 

In agreement with eq. (|14b|) the acting force F is also linear in its dependency on the 
radius Rob of the obstacle as drawn in fig. ^ (crosses and solid line). It can be expressed by 
the ffi 

F = (111 ■/?„, + 18) 10^^. (16) 

sec^ 

There exists a resulting force F ^ 18 ■ 10^ g ■ cm/sec for Rot -^ since as shown below 
even a very small obstacle hinders the stream particles (of finite radius) and results in a 
certain resistance for the stream. 

Thus, we find the surprising result that the functional form of the acting force as a 
function of the density of the stream p (eq. (|I^)) and of the radius Rob of the obstacle 
(eq. (|16])) which we measured in our simulations equals that of the crude approximation 
eq. (p^ . This coincidence is not trivial since the stream particles which collide with the 
obstacle build up a characteristic "corona" around the obstacle as it can be seen in fig. |^. 
The size and shape of this corona vary with the stream density p and the radius -Rob- 

For the case of a fiat obstacle one observes qualitatively the same behavior as for a 
spherical one. The snapshot fig. ^ shows that in this case the corona which forms the effective 
obstacle has approximately the same shape as for spherical obstacles. For the force as a 
function of the obstacle size we measure Ffl^^f ~ 156 ■ 10^ ■ Rob instead of F ~ 111 ■ 10^ ■ Rob for 
the simulation with a spherical obstacle of the same size and identical simulation parameters 
(fig. ^, dashed line). The shape of the obstacle causes only a change in the prefactor of the 
force law but leaves the functional form unaffected. 
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The color in figs. |^ and ^ codes for tlie absolute value of the particle's velocity. In both 
cases there is a wide region in front of the obstacle where the velocities are very small. In this 
region the particles almost agglomerate to the obstacle. The width of this region depends 
strongly on the stream velocity and on the material parameters. In figure |^ one observes 
a small gap of low particle density directly in front of the obstacle. In this gap one finds 
particles of high kinetic energy. The properties of the gap will be discussed below. 

The dependence of the force F on the velocity Vst of the stream is more complicate. 
Fig. P shows that the quadratic behavior from the estimation eq. (|14c|) is only valid for high 



velocities {vst > 700 cm/sec) in the case of a small obstacle. For lower velocities the curve 
deviates significantly from this function. The data can be fitted to 



,..1-5 



(17) 

however, we are not sure whether a power law is the adequate representation of these data. 
Hence, from eqs. (p^D, (|I6|) and (O) we find for our particular system 



. f ?y- Vst > 700 cm/sec 
p-Ro,-{ '' ' (18) 

vl:^ Vst < 700 cm/sec . 

For large impact velocity Vgt the corona is in direct contact with the obstacle as shown 
in the snapshots (figs. ^ and ^. In the region of the transition where the force starts to 
deviate from the quadratic law we observe that the agglomerated region of particles having 
low velocity in front of the obstacle becomes large. In figs. and |^ the size of this region is 
determined mainly by the size of the obstacle, whereas for small stream velocity the "cold" 
region may become significantly larger than the obstacle itself. Fig. ^ shows three snapshots 
of the system with stream velocity well above (left), at the transition point (middle) and 
well below the transition (right). 

Moreover close to the transition one observes a gap forming in between the corona and 



the obstacle. This gap is filled with granular gas at high temperature. In section |I^ we will 
discuss the field of granular temperature which supports this statement. 
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The width of this gap depends on the stream velocity as shown in fig. |TT]. For 
Vst ~ 700 cm/sec there is no gap. It appears at Vgt ^ 700 cm/sec and there is a sharp 
transition. The width of the gap increases for decreasing stream velocity. A similar effect, 
i.e. the appearance of a gap when varying the driving parameters of the system, was found 
recently in an event driven simulation of granular particles which are heated (with granular 
temperature) from below when gravity acts [|l^. Similar as in our case there seems to be a 
critical value for the driving parameter (in [0] the temperature) when the gap appears. 



For the region of parameters applied in our simulations the stream width b = 16cm does 
not influence the results. For obstacle size Rob = 0.5 cm, stream velocity Vst = 200 cm/sec 
and stream density p = 1 we made tests with bi = 2 b. We found for the time averaged 
force acting on the obstacle F{bi) = 1.009 F{b) which flts well into the region covered by 
the error bars. 

IV. CONTINUUM REPRESENTATION OF THE SCATTER EXPERIMENT 



It has been shown by several authors (e.g. [pr^-[T9|) that under certain conditions the 
dynamics of granular material can be described using hydrodynamics equations. In this 
section we want to present fleld data for the density p (r) and the temperature T (r) of the 
granular system. The fleld data have been generated by coarse graining the particle data 
to a grid with grid constant a = 0.1cm ^ R™°'^. The average has been taken over a time 
interval which corresponds to approximately 500 rows of stream particles colliding with the 
obstacle. 



Fig. |T2| shows the density and velocity flelds as color plots for Vst = 500cm/sec and three 
different radii of the spherical obstacle. Red color stands for high and blue for low values of 
the density or the velocity. All flgures in one row have the same scaling for the colors. For 
the case of the large sphere {Rob = 1.5cm) the density in front of the obstacle is very high. 
The particles which join this high density region are of low temperature. The temperature 
rises up to a characteristic distance from the surface of the obstacle and has its maximum 
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in the region where the incoming stream particles colhde with the corona which surrounds 
the obstacle, i.e. the location of maximum temperature coincides with the boundary of the 
corona. As discussed in the context of figs. and ^ for smaller obstacle size (or lower stream 
velocity) a region of low density and high temperature forms in front of the obstacle, i.e. the 
highest density can be found in a certain distance from the obstacle. The region of highest 
temperature does not necessarily coincide with the border of the corona. 

The fields of granular temperature and density as well as the appearance of the gap 
are reminiscent to the behavior of a bullet flying in air with supersonic velocity. It has 
been shown before that rapid granular flows may contain spatial regions in which the flow 
is supersonic ||T^j20| -p2| . The properties of the flow and the possible embedding into a 



continuum mechanics theory are subjects of work in progress [^]. 

V. DISCUSSION 

We investigated numerically the properties of a stream of granular particles which is scat- 
tered by an obstacle and the force acting upon the obstacle as a function of the stream density 
p, the stream velocity Vst and the obstacle size Rob using molecular dynamics technique. The 
fields of granular temperature and density have been generated from the molecular dynamics 
data by coarse graining. 

Surprisingly we find the force to be proportionally to the density p, the obstacle size Rob, 
and for high stream velocities proportionally to f^^, i.e. F ~ pRobvl^ as one would expect 
from a stream of non-interacting hard spheres. For lower velocities we find that the force 
deviates from the F ~ v^^ law, but still we find a linear law for the force as a function of 
Rob and p. This deviation occurs at a critical value for the stream velocity Vst- At the same 
velocity we observe a gap forming in between the surface of the obstacle and the corona of 
stream particles which agglomerate and rest around the obstacle. This gap rises when the 
stream velocity decreases. The widths of this gap which is filled with granular gas at high 
temperature grows with decreasing stream velocity. 
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This agglomeration of particles in front of the obstacle leads to an effective obstacle size 

eff 
RJ; larger than Rob- The shape of the corona does not sensitively depend on the shape of 

the obstacle, the functional form of the force F ~ p Rob v1^ remains invariant if the spherical 
obstacle is replaced by a fiat one. 

For the case of a sphere of radius i? in a viscous fluid one finds for very low stream 
velocity (low Reynolds number) experimentally the Stolies law F = Qthtj Rvgt- This re- 
sult can be derived analytically from the continuum mechanical Navier-Stokes equation for 
incompressible flow [^ 



p (vst ■ V) Vst = -Vp + riV^Vst (19) 

as result of the pressure in front and behind the obstacle. The conditions of low Reynolds 
number and high viscosity assure that streamlines of continuous media cling to each point 
of the obstacle, which is a condition for the derivation of the Stokes law. For the case of 
a granular material this condition does not hold, i.e. the pressure which acts from behind 
the obstacle is zero, or at least very close to zero since as visible from the snapshots figs. |^, 
M and |IU| and from the density fields in fig. O there are almost no particles which collide 



with the obstacle from behind. According to these considerations we conclude that one 
cannot expect Stokes law even for very low stream velocity Vst because the precondition of 
a homogeneous medium with high viscosity does not hold. 

The field data of the granular temperature and the density show that for large obstacle 
size the location of maximum temperature coincides with the boundary of the corona. For 
the case of smaller obstacles the highest density is found in a certain distance from the 
obstacle whereas the region close to the obstacle is filled by a diluted granular gas of high 
temperature. For fixed obstacle size this gap can be found below a certain stream velocity 
Vst- It grows rapidly for decreasing Vgt- 

Instead of using molecular dynamics technique for the simulation of the dynamics of the 
granular material which requires large computational effort, one could apply an "event driven 
algorithm" which works much faster than molecular dynamics. Using such algorithms one 

12 



does not solve Newton's equation of motion for each single collision between the grains but 
one calculates the relative tangential and normal velocities after each collision as functions 
of the relative impact velocities (normal and tangential). Therefore one needs the normal 
and tangential viscoelastic restitution coefficients which are functions of the material con- 
stants and the impact velocities. For the approximation that the dilatation is independent 
of the impact velocities Pao and others analytically calculated the normal restitution coef- 



ficient |2gj2^ . A more general approach for the derivation of the coefficients of restitution 



based on physical interaction forces between colliding spheres pj-Q has been given in P7 



Effective algorithms for these type of event driven molecular dynamics simulations can be 



found e.g. in |28-33 
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TABLES 



a = y/(l - 1/2) a = 10^° g • cm-°-5 • sec'^ 

j3 = A-a 15 = 1000 g • cm-^ • sec"^ 

shear friction 75 = 1000 sec~^ 

Coulomb friction parameter /x = 0.5 

material density p™ = 1 g/cm^ 

TABLE I. The parameters used in the simulations. 
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FIG. 1. The setup of the numerical experiment. 
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FIG. 2. The differential intensity of the stream scattered at a spherical obstacle over the angle 
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FIG. 3. The distribution of angular momentum L{Q) over the scatter angle G. 
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FIG. 4. The distribution of the absolute value of the averaged particle momentum over the 
scatter angle 0. Again, for a better statistics we took advantage of the symmetry of the problem. 
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FIG. 5. The force over the stream density p. 
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FIG. 6. The force over the obstacle size Rob- In agreement with eq. (|l4b| ) we find a linear 
function. The solid line shows the force experienced by a spherical obstacle, the dashed line shows 
the force for a flat obstacle (explanation see below). 
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FIG. 7. Snapshot of the simulation. The dissipation during the cohisions causes a stationary 
corona which shape and size depend on the parameters of the simulation. The color codes for the 
absolute value of the particle velocities, red means high velocity, blue means low. 
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FIG. 8. Snapshot of a simulation using a flat obstacle. Similar as in fig. 7 one finds a corona 
of stream particles shielding the obstacle. 
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FIG. 9. The force over the stream velocity Vgt- Small obstacle: for large velocity Vgt the force 
agrees with the crude approximation ( |14cD . For smaller impact velocity, however, the force deviates 
from this function. The dashed line shows the function F ~ v^''^. Large obstacle: the quadratic 
function remains valid even for smaller Vgt- We expect a deviation from this function for yet smaller 
impact velocity. 
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FIG. 10. For lower impact velocity u^t ^ gap in between the obstacle and the corona is observed. 
The gap appears at the same velocity when the force deviates from the quadratic law. Left: 
Vst = 2000 cm/sec, middle: Vst = 700 cm/sec, right: Vst = 200 cm/sec. Again the color codes 
for the absolute value of the particle velocities. The three figures are drawn with the same color 
scaling. 
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FIG. 11. The size of the gap as a function of the stream velocity. The value of Vgt when the 
gap appears coincides with the stream velocity when the force which acts on the obstacle deviates 



from the F ~ ii^^-behavior. 
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FIG. 12. The fields of density (top row) and granular temperature (bottom row) for three 
different obstacle sizes. Discussion see text. 
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